Visualización de datos geoespaciales con R.
En este taller veremos como visualizar datos geoespaciales (vectoriales) con R, sin embargo, los paquetes de R permiten hacer mucho más que sólo visualizar. Para aprender más sobre esto ver Geocomputation with R (Lovelace, Nowosad, y Muenchow 2019).
Instala los paquetes y descarga el material del taller.
install.packages(c("usethis", "tidyverse", "sf", "leaflet"))
usethis::use_course("https://github.com/tereom/taller-mapas-mv/blob/master/material.zip?raw=true")
Usaremos los paquetes tidyverse (Wickham et al. 2019) y sf(Pebesma 2018).
Objetivo: Un mapa de población por departamento de Uruguay
¿Qué necesito?
Los datos geoespaciales de los departamentos los obtenemos de mapas vectoriales del INE y usaremos la función st_read() para leerlos:
depart_uy <- st_read("datos/ine_depto")
Reading layer `ine_depto' from data source `/Users/tereom/Documents/mapas-taller-mv/datos/ine_depto' using driver `ESRI Shapefile'
Simple feature collection with 20 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 366582.2 ymin: 6127919 xmax: 858252.1 ymax: 6671738
Projected CRS: WGS 84 / UTM zone 21S
La función st_read() recibe como argumento la ruta a la carpeta (puede ser al archivo con terminación .shp) y la asigna a depart_uy.
Los podemos graficar directo con plot(depart_uy) o con ggplot (❤️) en el segundo caso usamos geom_sf() para agregar la capa de coordenadas:
ggplot(depart_uy) +
geom_sf()
Si inspeccionamos el objeto depart_uy descubrimos que es un poco familiar (y un poco nuevo), lo familiar es que es una tabla de datos de tipo data.frame, con una columna adicional de geometría que almacena la información de la geometría espacial.
class(depart_uy)
[1] "sf" "data.frame"
glimpse(depart_uy)
Rows: 20
Columns: 6
$ AREA_KM2_ <int> 530, 11928, 4536, 6106, 11643, 10417, 10016, 1392…
$ PERIMETER <dbl> 202136.96, 722385.94, 404534.80, 431056.43, 82132…
$ DEPTO <int> 1, 2, 3, 5, 6, 8, 9, 11, 12, 13, 14, 15, 16, 17, …
$ NOMBRE <chr> "MONTEVIDEO", "ARTIGAS", "CANELONES", "COLONIA", …
$ CDEPTO_ISO <chr> "UYMO", "UYAR", "UYCA", "UYCO", "UYDU", "UYFD", "…
$ geometry <MULTIPOLYGON [m]> MULTIPOLYGON (((579197 6136..., MULT…
En la gráfica podemos colorear (como en cuaqluier ggplot) los polígonos utilizando con el valor de alguna de las columnas del data.frame.
ggplot(depart_uy) +
geom_sf(aes(fill = AREA_KM2_))
También podemos unirle una tabla de datos externa, por ejemplo de población (datos de juntas departamentales) como lo haríamos de normal:
basicos <- read_csv("datos/datosbasicosjds.csv")
glimpse(basicos %>% select(1:10, 38))
Rows: 19
Columns: 11
$ `Junta Departamental` <chr> "Artigas", "Canelones", "Cerro Largo",…
$ `Tiene WEB` <chr> "NO", "Si", "Si", "Si", "Si", "Si", "S…
$ `Link de la web` <chr> NA, "http://www.juntadecanelones.gub.u…
$ Dirección <chr> "Cnel. Carlos Lecueder 510", "Luis Alb…
$ `Teléfono 1` <chr> "4772 4833", "4332 2420", "4642 2283",…
$ `Telefono 2` <chr> NA, NA, "4642 3471", NA, "4362 4778", …
$ Email <chr> "juntadeartigas@gmail.com", "juntacan@…
$ Plenarios <chr> "1er - 2do y 3er jueves de cada mes", …
$ `Hora Plenario` <chr> "20.30", "18.00", "19.00", "20.00", "1…
$ `Local propio` <chr> "Si", "Si", "Si", "Si", "Si", "No", "S…
$ Población <chr> "73.377", "520.173", "84.698", "123203…
# un poco de limpieza para poder unir las tablas
depart_pob <- basicos %>%
select(nombre = `Junta Departamental`, pob = `Población`) %>%
mutate(nombre = janitor::make_clean_names(nombre),
pob = as.numeric(str_remove_all(pob, "\\.")))
depart_uy_pob <- depart_uy %>%
mutate(nombre = janitor::make_clean_names(NOMBRE)) %>%
left_join(depart_pob)
glimpse(depart_uy_pob)
Rows: 20
Columns: 8
$ AREA_KM2_ <int> 530, 11928, 4536, 6106, 11643, 10417, 10016, 1392…
$ PERIMETER <dbl> 202136.96, 722385.94, 404534.80, 431056.43, 82132…
$ DEPTO <int> 1, 2, 3, 5, 6, 8, 9, 11, 12, 13, 14, 15, 16, 17, …
$ NOMBRE <chr> "MONTEVIDEO", "ARTIGAS", "CANELONES", "COLONIA", …
$ CDEPTO_ISO <chr> "UYMO", "UYAR", "UYCA", "UYCO", "UYDU", "UYFD", "…
$ nombre <chr> "montevideo", "artigas", "canelones", "colonia", …
$ pob <dbl> 1318755, 73377, 520173, 123203, 57084, 67047, 588…
$ geometry <MULTIPOLYGON [m]> MULTIPOLYGON (((579197 6136..., MULT…
Y graficamos nuestro mapa coloreado por población.
ggplot(depart_uy_pob) +
geom_sf(aes(fill = pob))
Ahora inspeccionamos el lado nuevo de nuestro objeto, cuando leímos los datos apuntamos a una carpeta:
depart_uy <- st_read("datos/ine_depto")
esta carpeta de donde leímos los datos tiene muchos archivos.
fs::dir_tree("datos/ine_depto/")
datos/ine_depto/
├── ine_depto.dbf
├── ine_depto.prj
├── ine_depto.sbn
├── ine_depto.sbx
├── ine_depto.shp
├── ine_depto.shx
└── ine_depto_metadatos.pdf
Esto es porque, en este ejemplo, nuestros datos vectoriales están almacenados en un shapefile (un formato muy común para datos geoespaciales).
Un shapefile es un grupo de archivos que contienen geometrías e información de acuerdo a un estándar especificado por el Insituto de Investigación en Sistemas de Ecosistemas (ESRI). En nuestro ejemplo tenemos los siguientes:
nombre.shp: es el archivo principal y contiene la geometría.
nombre.dbf: es la base de datos asoiciada al objeto y almacena la información de los atributos de los objetos, en el caso de ine_depto hay un renglón por polígono, es decir, cada línea contiene las características de un departamento.
AREA_KM2_ PERIMETER DEPTO NOMBRE CDEPTO_ISO
1 530 202137.0 1 MONTEVIDEO UYMO
2 11928 722385.9 2 ARTIGAS UYAR
3 4536 404534.8 3 CANELONES UYCA
4 6106 431056.4 5 COLONIA UYCO
5 11643 821325.4 6 DURAZNO UYDU
6 10417 669409.9 8 FLORIDA UYFD
crs_depart <- st_crs(depart_uy)
crs_depart$input
[1] "WGS 84 / UTM zone 21S"
nombre.shx: almacena el índice de las entidades geométricas.
nombre.sbn y .sbx (opcional): almacena índice espacial.
Como vimos en el ejemplo, read_sf() recibe la ruta a la carpeta o a alguno de los archivos de la capa de interés (en caso de tener varias capas en una carpeta podemos indicar la ruta hasta alguno de los archivos): read_sf("datos/ine_depto/ine_depto.shp") o utilizar el argumento layer: read_sf("datos/ine_depto", layer = "ine_depto")).
El segundo aspecto nuevo, es la columna geometry en los datos,
depart_uy$geometry
Geometry set for 20 features
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 366582.2 ymin: 6127919 xmax: 858252.1 ymax: 6671738
Projected CRS: WGS 84 / UTM zone 21S
First 5 geometries:
ésta puede almacenar distintos tipos de datos vectoriales.
Los tipos van aumentando en complejidad, el más simple serían puntos, cuya unión forma líneas, con estos dos formas polígonos. Algunas de las geometrías más comunes están representadas en la siguiente imagen:
Figure 1: Lovelace, Nowosad, y Muenchow (2019) (CC BY-NC-ND 4.0).
🗺 Grafica los polígonos de barrios de Montevideo (los archivos se encuentran en datos/ine_barrios).
🌎 Lee el archivo en datos/uptu_cultura y agrega las ubicaciones contenidas a la misma gráfica que los polígonos.
🌐 Describe la diferencia entre la geometría de los barrios de Montevideo y del archivo de cultura.
Notemos que podemos utilizar las operaciones usuales para transformar tablas de datos, por ejemplo, si queremos crear una columna de tipo de establecimiento usamos la función mutate de dplyr.
glimpse(cultura_mv)
Rows: 398
Columns: 3
$ NOMBRE <chr> "TEATRO ALFREDO MORENO", "TEATRO ALIANZA URUGUAY-E…
$ DIRECCION <chr> NA, "PARAGUAY 1217", "SAN JOSE 1426", "COLONIA 187…
$ geometry <POINT [m]> POINT (582125 6138279), POINT (573785.3 6136…
cultura_mv <- cultura_mv %>%
mutate(tipo_establecimiento = case_when(
str_detect(NOMBRE, "TEATRO") ~ "teatro",
str_detect(NOMBRE, "MUSEO") ~ "museo",
str_detect(NOMBRE, "BIBLIOTECA") ~ "biblioteca",
str_detect(NOMBRE, "CINE") ~ "cine",
TRUE ~ "otro"
))
glimpse(cultura_mv)
Rows: 398
Columns: 4
$ NOMBRE <chr> "TEATRO ALFREDO MORENO", "TEATRO ALIANZ…
$ DIRECCION <chr> NA, "PARAGUAY 1217", "SAN JOSE 1426", "…
$ geometry <POINT [m]> POINT (582125 6138279), POINT (57…
$ tipo_establecimiento <chr> "teatro", "teatro", "teatro", "teatro",…
Veamos que tenemos disponibles las opciones típicas de ggplot como facet_wrap().
ggplot() +
geom_sf(data = barrios_mv) +
geom_sf(data = cultura_mv, color = "red", alpha = 0.8, size = 0.5) +
labs(subtitle = "Establecimientos culturales en Montevideo",
color = "tipo") +
theme_void() +
facet_wrap(~tipo_establecimiento)
sf de una tabla de datosImaginemos que tenemos una tabla de datos georreferenciados, pueden ser datos colectados con gps en una encuesta o ubicaciones obtenidas de google maps. Lo podemos almacenar en un data.frame que podemos a su vez convertir en un objeto de tipo sf, lo que debemos hacer es indicar cuáles son las columnas que incluyen información geográfica.
library(googlesheets4)
url <- ("https://docs.google.com/spreadsheets/d/1h7n3dpjqGofWYtKf5ADCDAr8nT5U5UalMFq9nS3G_ew/edit#gid=0")
gs4_deauth()
colab <- read_sheet(url)
glimpse(colab)
Rows: 3
Columns: 3
$ nombre <chr> "Teresa", "María", "casa"
$ long <dbl> 18.93790, 19.36930, 18.95584
$ lat <dbl> -99.23250, -99.12646, -99.22654
Por ahora colab es simplemente un data.frame,
class(colab)
[1] "tbl_df" "tbl" "data.frame"
Lo podemos convertir a un objeto de clase sf con la función st_as_sf, indicando cuales son las columnas con información geográfica. De manera opcional podemos indicar también la información de la proyección (Google Maps utiliza la 4326 para reportar).
colab_longlat <- colab %>%
st_as_sf(coords = c("lat", "long")) %>%
st_sf(crs = 4326)
class(colab_longlat)
[1] "sf" "tbl_df" "tbl" "data.frame"
Y ahora lo podemos graficar con ggplot y geom_sf().
ggplot(colab_longlat) +
geom_sf(aes(color = nombre))
Como vimos con nuestros primeros mapas, los objetos espaciales suelen tener asociado un valor en sistema de referencia de coordenadas (CRS), en esta sección explicamos a qué se refiere y porque es importante.
La Tierra tiene forma elipsoidal y para ubicar puntos en la Tierra debemos especificar un modelo de elipse, un punto de origen y dirección de los ejes, y una estrategia para proyectar a un plano.
rgdal::projInfo(type = "ellps") %>%
rmarkdown::paged_table()
Los conceptos anteriores son la base de un sistema de coordenadas, es importante recalcar que no todos los datos geoespaciales usan el mismo sistema de coordenadas. Y para ser representados en un mismo mapa debemos cuadrarlos.
Figure 2: Comparación entre proyecciones tangente y secante cilíndirca, cónica y de azimut con licencia CC-BY 4.0
Las distorsiones resultan en pérdidas de información que pueden ser en área, forma, distancia o dirección. Por ejemplo Mercator preserva dirección y es útil para navegación:
states <- map_data("state")
usmap <- ggplot(states, aes(x=long, y=lat, group=group)) +
geom_polygon(fill="white", colour="black")
usmap + coord_map("mercator")
Por su parte Azimuthal Equal Area preserva area pero no dirección:
usmap + coord_map("azequalarea")
En particular en investigación es común usar Universal Transverse Mercator (UTM) porque preserva ángulos y dirección (por lo que es común en sistemas de navegación). Para minimizar la distorsión se divide la Tierra en 60 regiones y utiliza una proyección (de secantes) en cada una. Esta es la proyección más común para Uruguay.
Figure 3: Husos y Zonas UTM con licencia CC-BY 3.0
Como vimos arriba, los conjuntos de datos geoespaciales tienen asociado un sistema de referencia de coordenadas (y proyección si es el caso), Hay varias maneras de hacer referencia al sistema asociado a un conjunto de datos, una notación estándar que se utiliza para describir un CRS es proj4string y se ven como:
st_crs(barrios_mv)$proj4string # cartesianas
[1] "+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs"
st_crs(colab_longlat)$proj4string # geográficas
[1] "+proj=longlat +datum=WGS84 +no_defs"
También se puede usar el código ESPG asociado, los códigos EPSG. Los códigos EPSG (European Petroleum Survey Group) son códigos de 4 a 5 dígitos que representan definiciones de CRS:
st_crs(32721)$proj4string
[1] "+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs"
Hay más maneras de documentar el sistema de coordenadas, sin embargo las dos mencionadas son muy comunes.
🗾 En un editor de texto abre el archivo .prj de alguno de los conjuntos de datos geoespaciales que hemos utilizado.
🗾 ¿Cuáles hemos usado hasta ahora?
Muchas veces es necesario transformar nuestros datos a otro sistema de coordenadas geográficas, por ejemplo, si queremos trabajar con conjuntos de datos con distinto CRS, o para utilizar ciertas herramientas como leaflet.
Para hacer mapas interactivos usaremos el paquete leaflet (Cheng, Karambelkar, y Xie 2021):
Leaflet es una librería de JavaScript para realizar mapas interactivos, el paquete de R nos permite utilizarla desde R.
Leaflet recibe datos geoespaciales con CRS 4326 (al igual que las coordenadas que copiamos de GoogleMaps). Por ellos si queremos graficar los puntos culturales debemos comenzar reproyectando:
cultura_mv_longlat <- st_transform(cultura_mv, 4326)
Ahora graficamos, la sintaxis de leaflet es similar a la de ggplot, pero usamos %>% en lugar de +:
leaflet(data = cultura_mv_longlat) %>%
addTiles() %>%
addCircles()
La función addTiles() agrega las tejas del fondo. En este ejemplo agregagmos marcadores que muestran la ubicación de los sitios culturales de Montevideo, podemos añadir información de variables en los datos, por ejemplo, el nombre de la atracción. Y con addCircles() agregamos la geometría de puntos representados con círculos.
leaflet(data = cultura_mv_longlat) %>%
addTiles() %>%
addMarkers(popup = ~NOMBRE)
Con leaflet hay mucha flexibilidad, podemos meter html a las etiquetas o definir nuestros propios íconos, en el sitio se explica claramente con ejemplos.
rLadiesIcon <- makeIcon(
iconUrl = "https://raw.githubusercontent.com/rladies/starter-kit/master/logo/R-LadiesGlobal_CMYK_offline_LogoOnly.svg",
iconWidth = 38, iconHeight = 38
)
leaflet(data = colab_longlat) %>%
addTiles() %>%
addMarkers(popup = ~nombre, icon = rLadiesIcon)
También podemos agregar objetos espaciales de otros tipos, para polígonos se utiliza addPolygons().
barrios_longlat <- barrios_mv %>%
st_transform(crs = 4326)
leaflet(data = barrios_longlat) %>%
addTiles() %>%
addPolygons()
🌍 ¿Qué ocurre si utilizas leaflet con los datos no transformados, es decir si utilizas colab en lugar de colab_longlat?
🌍 Agrega al siguiente código la ubicación de los puntos de interés de Montevideo.
leaflet(data = colab_longlat) %>%
addTiles() %>%
addMarkers(popup = ~nombre, icon = rLadiesIcon)
🌍 Experimenta con otros fondos, utiliza la función addProviderTiles() y elige un fondo de aquí.
leaflet(data = colab) %>%
addProviderTiles(providers$Stamen.Watercolor) %>%
addMarkers(popup = ~nombre, icon = rLadiesIcon)
Las operaciones con datos geoespaciales van más allá del alcance de este taller, pero aquí hay algunos ejemplos de lo que se puede hacer.
barrios_mv %>%
mutate(area = st_area(geometry))
Simple feature collection with 62 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 551982.7 ymin: 6133499 xmax: 589227.2 ymax: 6159810
Projected CRS: WGS 84 / UTM zone 21S
# A tibble: 62 x 5
AREA_KM NOMBBARR NROBARRIO geometry area
* <dbl> <chr> <int> <MULTIPOLYGON [m]> [m^2]
1 2.11 Ciudad Vi… 1 (((573146.3 6137145, 573148.… 2297414…
2 1.30 Centro 2 (((574497.2 6137037, 574506.… 1294508…
3 0.694 Barrio Sur 3 (((574564 6136237, 574570.5 … 726727…
4 2.28 Cordon 4 (((574397.4 6137451, 574382.… 2278114…
5 0.798 Palermo 5 (((575149.1 6136881, 575144 … 832445…
6 0.766 Parque Ro… 6 (((576530.7 6135354, 576532.… 778529…
7 2.73 Punta Car… 7 (((576200.9 6135810, 576200.… 2818381…
8 4.16 Buceo 9 (((580161.9 6137234, 580167.… 4164074…
9 3.38 Pque. Bat… 10 (((578338.6 6137972, 578352.… 3371317…
10 3.50 Malvin 11 (((580642.4 6138757, 580590.… 3502102…
# … with 52 more rows
barrios_cultura_mv <- barrios_mv %>%
st_join(cultura_mv, join = st_intersects) # usualmente st_within es más apropiado
barrios_cultura_mv %>%
st_drop_geometry() %>%
group_by(NOMBBARR) %>%
summarise(n_lugares = sum(!is.na(NOMBRE))) %>%
arrange(-n_lugares)
# A tibble: 62 x 2
NOMBBARR n_lugares
<chr> <int>
1 Centro 81
2 Cordon 57
3 Ciudad Vieja 50
4 Pque. Batlle, V. Dolores 24
5 Buceo 14
6 Palermo 13
7 Punta Carretas 13
8 Parque Rodo 10
9 Pocitos 10
10 Aguada 8
# … with 52 more rows
colab_longlat %>%
st_distance()
Units: [m]
[,1] [,2] [,3]
[1,] 0.000 49039.48 2082.469
[2,] 49039.476 0.00 46963.494
[3,] 2082.469 46963.49 0.000
La cheatsheet del paquete sf describe operaciones geométricas disponibles, e incluye útiles esquemas.
Texto y gráficas con licencia Creative Commons Attribution CC BY 4.0, a menos que se indique lo contrario o tengan otra fuente (indicado en la nota). Código en https://github.com/tereom/taller-mapas-mv.